1 Introducció

Aquest és un exercici de mineria de dades sobre l’estudi dels terratrèmols o activitat sísmica global, és a dir plantegem l’estudi de l’activitat sísmica mundial.

Disposar d’un bon conjunt de dades i un bon coneixement del camp del que volem fer l’estudi ajuda a fer-se bones preguntes i a orientar-se cap a l’objectiu on es vol arribar. En aquest cas, no disposo d’un bon coneixement de geologia i les dades…ja veurem, però sí que personalment, crec és un tema interessant i val la pena dedicar-li aquesta pràctica.

2 Elecció del conjunt de dades

Comencem doncs amb un joc de dades procedents de Kaggle, en particular:

Aquest dataset que hem triat, Global Significant Earthquake Database from 2150BC, és un llistat d’uns 6000 terratrèmols que van del 2150 BC fins a l’actualitat. Això ja ens diu de bones a primeres que les dades prèvies a l’ús d’instruments i recull de dades, probablement són incompletes. De fet, ja avanço que aquest conjunt de dades és altament incomplet.

3 Objectius

Com a objectiu general, volem mostrar quines són les zones i països amb més activitat sísmica i perill de tsunamis i el seu impacte en aquestes àrees, això posarà de manifest l’estructura, posició i dinàmica de les plaques tectòniques i a la vegada les possibles conseqüències per la població. Volem veure quins països són els més actius sísmicament. També volem fer una classificació i veure què fa que un terratrèmol sigui devastador, (significant). Intentarem esbrinar relacions o correlacions entre magnitud del sisme, profunditat a la que es produeix i altres variables. El valor d’aquestes dades pot ajudar a minimitzar conseqüències o a preveure i millorar infraestructures. De totes maneres com s’ha vist sempre, les àrees amb menys recursos i més pobres són les més damnificades.

4 Inspecció de les dades

library(tidyverse)

El primer que fem és carregar les dades i fer una ullada per veure de quantes variable, de quin tipus són i de quants registres disposem.

dades <- read_csv('Worldwide-Earthquake-database.csv')
## Rows: 6193 Columns: 47
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (4): FLAG_TSUNAMI, COUNTRY, STATE, LOCATION_NAME
## dbl (43): I_D, YEAR, MONTH, DAY, HOUR, MINUTE, SECOND, FOCAL_DEPTH, EQ_PRIMA...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
glimpse(dades)
## Rows: 6,193
## Columns: 47
## $ I_D                                <dbl> 1, 2, 3, 5877, 8, 11, 9712, 12, 13,…
## $ FLAG_TSUNAMI                       <chr> "No", "Yes", "No", "Yes", "No", "No…
## $ YEAR                               <dbl> -2150, -2000, -2000, -1610, -1566, …
## $ MONTH                              <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ DAY                                <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ HOUR                               <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ MINUTE                             <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ SECOND                             <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ FOCAL_DEPTH                        <dbl> NA, NA, 18, NA, NA, NA, NA, NA, NA,…
## $ EQ_PRIMARY                         <dbl> 7.3, NA, 7.1, NA, NA, NA, NA, 6.5, …
## $ EQ_MAG_MW                          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ EQ_MAG_MS                          <dbl> NA, NA, 7.1, NA, NA, NA, NA, NA, NA…
## $ EQ_MAG_MB                          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ EQ_MAG_ML                          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ EQ_MAG_MFA                         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ EQ_MAG_UNK                         <dbl> 7.3, NA, NA, NA, NA, NA, NA, 6.5, 6…
## $ INTENSITY                          <dbl> NA, 10, 10, NA, 10, 10, NA, NA, NA,…
## $ COUNTRY                            <chr> "JORDAN", "SYRIA", "TURKMENISTAN", …
## $ STATE                              <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ LOCATION_NAME                      <chr> "JORDAN:  BAB-A-DARAA,AL-KARAK", "S…
## $ LATITUDE                           <dbl> 31.100, 35.683, 38.000, 36.400, 31.…
## $ LONGITUDE                          <dbl> 35.50, 35.80, 58.20, 25.40, 35.30, …
## $ REGION_CODE                        <dbl> 140, 130, 40, 130, 140, 130, 140, 1…
## $ DEATHS                             <dbl> NA, NA, 1, NA, NA, NA, NA, NA, NA, …
## $ DEATHS_DESCRIPTION                 <dbl> NA, 3, 1, NA, NA, NA, NA, NA, NA, N…
## $ MISSING                            <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ MISSING_DESCRIPTION                <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ INJURIES                           <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ INJURIES_DESCRIPTION               <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ DAMAGE_MILLIONS_DOLLARS            <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ DAMAGE_DESCRIPTION                 <dbl> 3, NA, 1, NA, 3, NA, 3, 3, 3, 3, NA…
## $ HOUSES_DESTROYED                   <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ HOUSES_DESTROYED_DESCRIPTION       <dbl> NA, NA, 1, NA, NA, NA, NA, NA, NA, …
## $ HOUSES_DAMAGED                     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ HOUSES_DAMAGED_DESCRIPTION         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TOTAL_DEATHS                       <dbl> NA, NA, 1, NA, NA, NA, NA, NA, NA, …
## $ TOTAL_DEATHS_DESCRIPTION           <dbl> NA, 3, 1, 3, NA, NA, NA, NA, NA, NA…
## $ TOTAL_MISSING                      <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TOTAL_MISSING_DESCRIPTION          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TOTAL_INJURIES                     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TOTAL_INJURIES_DESCRIPTION         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TOTAL_DAMAGE_MILLIONS_DOLLARS      <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TOTAL_DAMAGE_DESCRIPTION           <dbl> NA, NA, 1, 3, NA, NA, 3, NA, NA, NA…
## $ TOTAL_HOUSES_DESTROYED             <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TOTAL_HOUSES_DESTROYED_DESCRIPTION <dbl> NA, NA, 1, NA, NA, NA, NA, NA, NA, …
## $ TOTAL_HOUSES_DAMAGED               <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TOTAL_HOUSES_DAMAGED_DESCRIPTION   <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…

Com podem veure, tenim 6193 registres i 47 variables. Aquestes dades, inclouen dades que es remonten a l’any 2150 BC fins a l’actualitat. Observem la quantitat de NA’s que hi ha!.

Fem una descripció d’aquestes variables, (https://medium.com/@dmukherjeetextiles/significant-earthquakes-2150-bc-to-2023-ad-e550c60fe863):

DIMENSIONS TEMPORALS

DIMENSIONS QUE DESCRIUEN LES CARACTERÍSTIQUES

DIMENSIONS GEOGRÀFIQUES

DIMENSIONS QUE INDIQUEN CONSEQÜÈNCIES HUMANES I MATERIALS

La resta són totals que de valors que ja s’han descrit anteriorment.

Les dades que farem servir un moment o altre són: YEAR, MONTH, DAY, HOUR, FLAG_TSUNAMI, FOCAL_LENGTH, EQ_PRIMARY, INTENSITY, COUNTRY, LOCATION_NAME, LATITUDE, LONGITUDE, DEATHS, DEATHS_DESCRIPTION, DAMAGE_MILLIONS_DOLLARS, DAMAGE_DESCRIPTION.

5 Preprocessament i gestió de característiques

Podem fer una primera ullada i mirar com es distribueixen les dades recollides al llarg dels anys.

ggplot(dades, 
       aes(x = YEAR)
       ) + 
  geom_histogram(binwidth = 2,
           color = "darkred",
           fill = "darkred" )

Podem veure que en temps més recents tenim més dades. Això és raonable ja que a mida que tirem enrere tenim menys informació de les variables.

5.1 Neteja

Si mirem els valors NA:

print('NA')
## [1] "NA"
colSums(is.na(dades))
##                                I_D                       FLAG_TSUNAMI 
##                                  0                                  0 
##                               YEAR                              MONTH 
##                                  0                                407 
##                                DAY                               HOUR 
##                                561                               2042 
##                             MINUTE                             SECOND 
##                               2247                               3364 
##                        FOCAL_DEPTH                         EQ_PRIMARY 
##                               2965                               1791 
##                          EQ_MAG_MW                          EQ_MAG_MS 
##                               4872                               3265 
##                          EQ_MAG_MB                          EQ_MAG_ML 
##                               4391                               6009 
##                         EQ_MAG_MFA                         EQ_MAG_UNK 
##                               6179                               5416 
##                          INTENSITY                            COUNTRY 
##                               3378                                  0 
##                              STATE                      LOCATION_NAME 
##                               5872                                  1 
##                           LATITUDE                          LONGITUDE 
##                                 54                                 50 
##                        REGION_CODE                             DEATHS 
##                                  1                               4129 
##                 DEATHS_DESCRIPTION                            MISSING 
##                               3649                               6172 
##                MISSING_DESCRIPTION                           INJURIES 
##                               6172                               4955 
##               INJURIES_DESCRIPTION            DAMAGE_MILLIONS_DOLLARS 
##                               4769                               5685 
##                 DAMAGE_DESCRIPTION                   HOUSES_DESTROYED 
##                               1762                               5411 
##       HOUSES_DESTROYED_DESCRIPTION                     HOUSES_DAMAGED 
##                               4495                               5710 
##         HOUSES_DAMAGED_DESCRIPTION                       TOTAL_DEATHS 
##                               5265                               4500 
##           TOTAL_DEATHS_DESCRIPTION                      TOTAL_MISSING 
##                               4164                               6168 
##          TOTAL_MISSING_DESCRIPTION                     TOTAL_INJURIES 
##                               6167                               4940 
##         TOTAL_INJURIES_DESCRIPTION      TOTAL_DAMAGE_MILLIONS_DOLLARS 
##                               4760                               5740 
##           TOTAL_DAMAGE_DESCRIPTION             TOTAL_HOUSES_DESTROYED 
##                               2920                               5380 
## TOTAL_HOUSES_DESTROYED_DESCRIPTION               TOTAL_HOUSES_DAMAGED 
##                               4416                               5772 
##   TOTAL_HOUSES_DAMAGED_DESCRIPTION 
##                               5384

Podem veure que hi ha moltíssims valors NA, algunes variables d’entrada són inservibles, d’altres pràcticament el 50% dels seus valors són nuls, d’altres ronda el 20% i poques es podríen netejar directament. Només FLAG_TSUNAMI, YEAR i COUNTRY estan netes, (3 de 47 !). Amb aquest panorama, si féssim una neteja com cal ens quedaríem a les fosques. Per tant decidim netejar algunes variables i conviure amb NA’s d’algunes altres, ho veurem tot seguit.

Abans que res, eliminem les darreres 16 variables que no farem servir per res i estan plegades de valors nuls.

dades <- dades |> 
  select(1:31)

Després d’anar comprovant com disminueix el nombre de registres(files), fem el següent:

  • Eliminem les dades anteriors al 1900. És a dir, considerem l’activitat sísmica del 1900 endavant.
  • Tenint en compte les variables que més ens interessen, netegem els camps: FOCAL_DEPTH (profunditat), EQ_PRIMARY (magnitud), HOUR (hora) i DEATHS (nombre de morts).
dades_1900 <- dades |> 
  filter(YEAR >= 1900)

dades_netes <- dades_1900 |>
  filter(YEAR >= 1900 & !is.na(FOCAL_DEPTH) & !is.na(EQ_PRIMARY) & !is.na(HOUR) & !is.na(DEATHS) & !is.na(DAMAGE_DESCRIPTION))

Podem fer un resum de com ens ha quedat el nostre dataset, (ho farem usant la llibreria skimr,
https://search.r-project.org/CRAN/refmans/skimr/html/skim.html).

library(skimr)

skim_without_charts(dades_netes)
Data summary
Name dades_netes
Number of rows 1135
Number of columns 31
_______________________
Column type frequency:
character 4
numeric 27
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
FLAG_TSUNAMI 0 1.00 2 3 0 2 0
COUNTRY 0 1.00 3 36 0 99 0
STATE 1108 0.02 2 2 0 8 0
LOCATION_NAME 0 1.00 4 52 0 994 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100
I_D 0 1.00 6096.53 2291.86 2583.00 4569.00 5291.00 7674.50 10494.00
YEAR 0 1.00 1984.63 27.54 1900.00 1969.50 1990.00 2006.00 2020.00
MONTH 0 1.00 6.55 3.42 1.00 4.00 7.00 9.00 12.00
DAY 0 1.00 16.09 8.59 1.00 9.00 16.00 24.00 31.00
HOUR 0 1.00 11.32 7.11 0.00 5.00 11.00 18.00 23.00
MINUTE 3 1.00 28.92 17.42 0.00 14.00 28.00 43.00 59.00
SECOND 113 0.90 30.49 17.31 0.10 15.00 30.20 45.80 59.90
FOCAL_DEPTH 0 1.00 31.67 40.55 0.00 10.00 22.00 33.00 651.00
EQ_PRIMARY 0 1.00 6.38 0.97 2.10 5.70 6.40 7.10 9.50
EQ_MAG_MW 597 0.47 6.54 0.88 4.20 5.90 6.40 7.18 9.50
EQ_MAG_MS 356 0.69 6.48 0.97 3.50 5.80 6.50 7.30 8.80
EQ_MAG_MB 422 0.63 5.90 0.64 2.10 5.50 5.90 6.30 8.20
EQ_MAG_ML 1100 0.03 5.96 1.07 3.60 5.45 6.00 6.65 7.70
EQ_MAG_MFA 1134 0.00 6.30 NA 6.30 6.30 6.30 6.30 6.30
EQ_MAG_UNK 983 0.13 6.34 0.87 3.60 5.80 6.30 7.03 8.30
INTENSITY 636 0.44 7.89 1.85 2.00 7.00 8.00 9.00 12.00
LATITUDE 0 1.00 19.68 21.46 -54.00 2.59 27.38 37.15 61.02
LONGITUDE 0 1.00 38.36 79.62 -178.25 -0.97 52.19 103.37 178.29
REGION_CODE 0 1.00 102.80 54.92 10.00 40.00 130.00 150.00 170.00
DEATHS 0 1.00 1875.49 15372.49 1.00 2.00 8.00 60.00 316000.00
DEATHS_DESCRIPTION 0 1.00 1.56 1.01 1.00 1.00 1.00 2.00 4.00
MISSING 1116 0.02 2410.89 9946.70 1.00 6.50 21.00 157.00 43476.00
MISSING_DESCRIPTION 1116 0.02 1.84 1.01 1.00 1.00 1.00 3.00 4.00
INJURIES 422 0.63 3676.10 34627.61 1.00 21.00 100.00 400.00 799000.00
INJURIES_DESCRIPTION 339 0.70 2.26 1.13 1.00 1.00 2.00 3.00 4.00
DAMAGE_MILLIONS_DOLLARS 758 0.33 1652.10 7791.53 0.10 5.00 35.00 400.00 100000.00
DAMAGE_DESCRIPTION 0 1.00 2.51 1.05 1.00 2.00 2.00 3.00 4.00

Les dades s’han reduït de 6193 a 1135 registres, (això és considerable i potser segons com pot donar lloc a dades esbiaxades, però pels nostres objectius considerem que no).

Nota: les dades que tenim del 1900 endavant (dades_1900), sense netejar són 3720 registres, les farem servir en algunes visualitzacions, així tindrem més punts representatius, ara bé pel cas de l’estudi PCA necessitarem les dades amb NA’s a zero ja que sinó tenim problemes.

Per les dades que es mostren no s’observen outliers, en el sentit de dades fora del seu rang, i tampoc tenim dades repetides. Podem visualitzar les distribucions de diferents variables:

library(cowplot)
anys <- ggplot(dades_netes, 
                    aes(x = YEAR)) +
  geom_bar(color = "black", fill = "cornflowerblue")

mesos <- ggplot(dades_netes, aes(x = MONTH)) +
  geom_bar(color = "black", fill = "cornflowerblue")

plot_grid(anys, mesos)

dies <- ggplot(dades_netes, 
                    aes(x = DAY)) +
  geom_bar(color = "black", fill = "cornflowerblue")

hores <- ggplot(dades_netes, aes(x = HOUR)) +
  geom_bar(color = "black", fill = "cornflowerblue")

plot_grid(dies, hores)

profunditat <- ggplot(dades_netes, 
                    aes(x = FOCAL_DEPTH)) +
  geom_bar(color = "black", fill = "cornflowerblue")

magnitud <- ggplot(dades_netes, aes(x = EQ_PRIMARY)) +
  geom_bar(color = "black", fill = "cornflowerblue")

plot_grid(profunditat, magnitud)

morts <- ggplot(dades_netes, 
                    aes(x = DEATHS)) +
  geom_bar(color = "black", fill = "cornflowerblue")

danys <- ggplot(dades_netes, aes(x = DAMAGE_MILLIONS_DOLLARS)) +
  geom_bar(color = "black", fill = "cornflowerblue")

plot_grid(morts, danys)

Podem veure com més cap a l’actualitat tenim més dades, no s’observa cap mes, dia o hora preferents en que es produeixin més terratrèmols, (es podria mirar de veure si el nombre de morts estan relacionats amb l’hora del dia en que s’ha produït el sisme).
En quant a la profunditat i magnitud, podria semblar que tenim outliers, però els valors estan a dins dels establerts i a part recordem que hem fet una neteja exhaustiva.

5.2 Anàlisi i visualització

La distribució de terratrèmols, tenint en compte la magnitud. (Ho representem en un mapa, per fer-lo més gran fem servir:
https://ggplot2-book.org/maps.html, https://www.tidyverse.org/blog/2020/08/taking-control-of-plot-scaling/) :

library(maps)
world <- map_data('world') 

title <- "Mapa de Terratrèmols segons la magnitud"
wr <- map_data("world")

p <- ggplot(dades_1900) + 
  geom_map(data = world, map = world, 
           aes(map_id = region), 
           fill="white", 
           colour="#7f7f7f") +
  expand_limits(x = world$long, y = world$lat) +
  geom_point(aes(x = LONGITUDE, y = LATITUDE, color = EQ_PRIMARY),
             size = 1) + 
  ggtitle(title) +
  labs(x = "Longitud", y = "Latitud", color = "Magnitud") +
  scale_colour_gradient(high = "#132B43", low = "lightblue")

pngfile <- fs::path(knitr::fig_path(),  "mapa1.png")

library(ragg)
agg_png(pngfile, width = 25, height = 12, units = "cm", res = 300)
plot(p)
## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
invisible(dev.off())
knitr::include_graphics(pngfile)

Aquest mapa mostra l’activitat sísmica més significativa. Es pot veure com els punts estan alineats indicant els límits de les plaques tectòniques.

Considerem ara la presència de tsunamis:

dades_1900$tsunami <- ifelse(dades_1900$FLAG_TSUNAMI %in% c("Yes"), "Yes", "No")
counts <- table(dades_1900$tsunami)

barplot(prop.table(counts), 
        col = c("cornflowerblue","darkred"), 
        main="Terratèmols i tsunamis", 
        legend.text=c("No Tsunami","Sí Tsunami"),
        xlab ="Tsunamis", ylab = "Percentatge",
        ylim=c(0,0.8))

S’observa que un 75% dels terratrèmols de les nostres dades no provoquen tsunamis, mentre que un 25% sí. Això pot indicar que una gran part dels sismes es produeixen a zones terrestres i l’altre restant en zones costeres o marítimes.

Això ho podem representar de nou en un mapa.

world <- map_data('world') 

title <- "Mapa de Terratrèmols segons la presència de tsunamis"
wr <- map_data("world")

p <- ggplot(dades_1900) + 
  geom_map(data = world, map = world, 
           aes(map_id = region), 
           fill="white", 
           color="#7f7f7f") +
  expand_limits(x = world$long, y = world$lat) +
  geom_point(aes(x = LONGITUDE, y = LATITUDE, color = FLAG_TSUNAMI),
             size = 1) + 
  ggtitle(title) +
  labs(x = "Longitud", y = "Latitud", color = "Hi ha hagut tsunami?")

pngfile <- fs::path(knitr::fig_path(),  "mapa2.png")

agg_png(pngfile, width = 25, height = 12, units = "cm", res = 300)
plot(p)
## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
invisible(dev.off())
knitr::include_graphics(pngfile)

Podem mirar les morts provocades pels terratrèmols que produeixen tsunami o no.

death_ts <- dades_1900 |> 
  filter(FLAG_TSUNAMI == "Yes")
death_nts <- dades_1900 |> 
  filter(FLAG_TSUNAMI == "No")
par(mfcol = c(1, 2))
hist(death_ts$DEATHS, breaks = 50, main = "Morts amb tsunami", xlab = "Nombre de morts", col = "cornflowerblue")
hist(death_nts$DEATHS, breaks = 50, main = "Morts sense tsunami", xlab = "Nombre de morts", col = "cornflowerblue")

Com és d’esperar, segons els resultats anteriors, les morts per sismes sense tsunami són d’un ordre de magnitud més gran. S’observa també que tot i que els grans terratrèmols provoquen gran quantitat de morts, (tant si hi ha tsunami o no), sempre hi ha un mínim de morts acumulats a quantitats més petites.

Si mirem els morts considerant la magnitud del terratrèmol:

graph1 <- ggplot(dades_netes,
       aes(x = EQ_PRIMARY, y = DEATHS)) +
  geom_point(color = "darkred") + 
  labs(x = "Magnitud", y = "Nombre de morts")

graph2 <- ggplot(dades_netes,
       aes(x = EQ_PRIMARY, y = DEATHS_DESCRIPTION)) +
  geom_point(color = "darkred") + 
  labs(x = "Magnitud", y = "Grup de mortalitat")
# https://ggplot2-book.org/arranging-plots
library(patchwork)
graph1 + graph2

Podem veure, complementant el resultat anterior, que la mortalitat s’incrementa amb la magnitud del terratrèmol, tot i que aquest resultat també posa de manifest que la densitat de població és un factor que incideix en la mortalitat, això es pot veure amb la presència de sismes de gran magnitud i pocs morts. (Recordem que els grups de mortalitat són: 0 = Cap, 1 = Poques, 2 = Algunes, 3 = Moltes i 4 = Moltíssims).

Podem mirar la magnitud dels terratrèmols que provoquen tsunami i els que no en provoquen.

mag_tsu <- dades_netes  |>
 filter(FLAG_TSUNAMI == "Yes")

mag_no_tsu <- dades_netes  |>
 filter(FLAG_TSUNAMI == "No")

g1 <- ggplot(mag_tsu, 
        aes( x = EQ_PRIMARY)) + 
   geom_bar(color = "black", fill = "cornflowerblue") +
   labs(title = "Magnitud amb tsunami", 
        x = "Magnitud")

g2 <- ggplot(mag_no_tsu, 
        aes( x = EQ_PRIMARY)) + 
   geom_bar(color = "black", fill = "cornflowerblue") +
   labs(title = "Magnitud sense tsunami", 
        x = "Magnitud")

plot_grid(g1, g2)

Es veu que els terratrèmols que provoquen tsunami tenen magnituds cap a la banda alta, els que no provoquen tsunami hi tenim un ventall més ampli de magnituds, (tenint en compte amb les dades amb les que treballem).

Centrem-nos ara amb la relació entre magnitud i profunditat:

ggplot(dades_1900,
       aes(x = FOCAL_DEPTH, y = EQ_PRIMARY, color = EQ_PRIMARY)) +
  geom_point() +
  labs(x = "Profunditat (km)", y = "Magnitud", color = "Magnitud")

S’observa la presència de terratrèmols de totes les magnituds per profunditats baixes (< 100 km), aquests són la majoria , però a partir dels 200 km de profunditat les magnituds són superiors a 6, i a partir dels 300 km la magnitud és superior a 7. (Caldria fer una mica de recerca per trobar alguna relació amb l’estructura de capes terrestres i límits de plaques tectòniques).

Considerem els terratrèmols més profunds (profunditat > 200 km) i situem-los al mapa.

earthq300 <- dades_1900 |> 
  filter(FOCAL_DEPTH >= 200)

world <- map_data('world') 

title <- "Terratrèmols més profunds ( > 200 km)"
wr <- map_data("world")

p <- ggplot(earthq300) + 
  geom_map(data = world, map = world, 
           aes(map_id = region), 
           fill="white", 
           color="#7f7f7f") +
  expand_limits(x = world$long, y = world$lat) +
  geom_point(aes(x = LONGITUDE, y = LATITUDE, color = EQ_PRIMARY),
             size = 3) + 
  ggtitle(title) +
  labs(x = "Longitud", y = "Latitud", color = "Magnitud") +
  scale_colour_gradient(high = "#132B43", low = "lightblue")

pngfile <- fs::path(knitr::fig_path(),  "mapa4.png")

agg_png(pngfile, width = 25, height = 12, units = "cm", res = 300)
plot(p)
invisible(dev.off())
knitr::include_graphics(pngfile)

Hem anat fent visualitzacions en el mapa, però ara volem especificar la relació dels terratrèmols amb la geografia, és a dir països on n’hi ha més i a on són més greus.

Comencem amb amb els països on l’activitat sísmica és més gran, llistarem els 20 primers:

top20_quant <- dades_1900  |>
  group_by(COUNTRY) |> 
  summarize(earthq_num = n()) |> 
  filter(COUNTRY != "")  |> 
  arrange(desc(earthq_num)) |> 
  head(n = 20)

# https://sebastiansauer.github.io/ordering-bars/
ggplot(top20_quant, 
       aes(x = reorder(COUNTRY, -earthq_num), y = earthq_num, 
           fill = earthq_num)) + 
  geom_bar(position = "dodge", stat = "identity") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  labs(title = "Països amb més nombre de terratrèmols", 
       x = "País", y = "Nombre de terratrèmols",
       fill = "Nombre terratrèmols") +
  scale_fill_gradient(high = "#132B43", low = "lightblue")

Volem ara llistar els 20 pitjors terratrèmols en magnitud, a quines zones han passat:

top20_mag <- dades_1900  |>
 group_by(LOCATION_NAME , EQ_PRIMARY) |> 
 filter(LOCATION_NAME != "") |> 
 arrange(desc(EQ_PRIMARY)) |> 
 head(n = 20)

 ggplot(top20_mag, 
        aes( x = reorder(LOCATION_NAME, -EQ_PRIMARY), y = EQ_PRIMARY,
             fill = EQ_PRIMARY)) + 
   geom_bar(position = "dodge", stat = "identity") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  labs(title = "Terratrèmols més forts", 
       x = "Zona", y = "Magnitud",
       fill = "Magnitud") +
  scale_fill_gradient(high = "#132B43", low = "lightblue")

Els terratrèmols que han causat més morts, (que no han de coincidir amb els anteriors, és a dir amb els de més magnitud).

top20_dead <- dades_1900  |>
 group_by(LOCATION_NAME , DEATHS) |> 
 filter(LOCATION_NAME != "") |> 
 arrange(desc(DEATHS)) |> 
 head(n = 20)

 ggplot(top20_dead, 
        aes( x = reorder(LOCATION_NAME, -DEATHS), y = DEATHS, fill = DEATHS)) + 
   geom_bar(position = "dodge", stat = "identity") +
   theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
   labs(title = "Terratrèmols més mortals", 
        x = "Zona", y = "Morts",
        fill = "Morts") + 
   scale_fill_gradient(high = "#132B43", low = "lightblue")

I els més destructius, (danys materials):

top20_dam <- dades_1900  |>
 group_by(LOCATION_NAME , DAMAGE_DESCRIPTION) |> 
 filter(LOCATION_NAME != "") |> 
 arrange(desc(DAMAGE_DESCRIPTION)) |> 
 head(n = 20)

 ggplot(top20_dam, 
        aes( x = reorder(LOCATION_NAME, -DAMAGE_DESCRIPTION), y = DAMAGE_DESCRIPTION, fill = DAMAGE_DESCRIPTION)) + 
   geom_bar(position = "dodge", stat = "identity") +
   theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
   labs(title = "Terratrèmols més detructius", 
        x = "Zona", y = "Index de destrucció",
        fill = "Index de destrucció") + 
   scale_fill_gradient(high = "#132B43", low = "lightblue")

Observacions:

D’aquestes darreres 4 gràfiques es poden observar algunes característiques. Xina encapçala els països amb més nombre de terratrèmols, Japó també. En quant a magnitud, destaquem Xile, Indonèsia, Japó, Rússia (Kamtxatka) i Alaska. Si mirem el nombre de morts, Haití és el primer de la llista i no apareix a les llistes anteriors, això ens pot indicar que és una zona altament poblada i un país molt pobre. En aquesta llista de més mortalitat, no hi apareixen ni Kamtxatka ni Alaska (amb terratrèmols de gran magnitud), ja que són zones amb una densitat de població molt baixa.
En quant al llistat de la destrucció, llevat del terratrèmol de San Francisco a California a principis del segle passat, s’observa que en general els països que hi apareixen són majoritàriament països amb construccions més precàries, destacant però Japó que està més enrere, que com és conegut, és un dels països amb les millors mesures antisísmiques.

Podem completar aquestes visualitzacions la relació entre el grau de danys relacionat amb la magnitud.

# https://r-graph-gallery.com/2d-density-plot-with-ggplot2.html
# https://stackoverflow.com/questions/43515112/reversing-default-scale-gradient-ggplot2
top20_dam <- dades_1900  |>
  group_by(DAMAGE_DESCRIPTION)

ggplot(dades_1900, aes(x=EQ_PRIMARY, y=DAMAGE_DESCRIPTION)) +
  geom_bin2d(bins = 10) +
  scale_fill_continuous(high = "royalblue4", low = "steelblue1") +
  labs(x = "Magnitud" , y = "Grau de danys")
## Warning: Removed 1145 rows containing non-finite outside the scale range
## (`stat_bin2d()`).

D’aquesta visualització no podem treure gaires conclusions, ja hem vist que els danys i la magnitud depèn molt de la densitat de població i de la riquesa dels països, (primer món vs tercer món).

5.3 Codificació

Creem una nova variable numèrica a partir de FLAG_TSUNAMI, que pendrà valor 1 o 0 segons sigui “Yes” o “No”

# Convertim Yes = 1, No = 0
dades_netes <- dades_netes |> 
  mutate(tsunami = ifelse(FLAG_TSUNAMI %in% c("No"), 0, 1))

5.4 Correlacions

Tot seguit cercarem correlacions en funció dels morts i altres variables que poden estar relacionades, (inclosa tsunami que acabem de crear).

# Utilitzem aquesta llibreria per fer servir la funcio multiplot()
library('Rmisc')
#Crearem primer una llista per mostrar les gràfiques de correlacions
#Crearem una llista per mostrar els atributs que interessen.

n = c("YEAR", "DAY", "HOUR", "tsunami", "FOCAL_DEPTH", "INTENSITY", "LATITUDE", "LONGITUDE", "EQ_PRIMARY", "DAMAGE_DESCRIPTION", "REGION_CODE") 
dataAux <- dades_netes |> 
  select(all_of(n))
histList2<- vector('list', ncol(dataAux))
for(i in seq_along(dataAux)){
  message(i)
histList2[[i]]<-local({
  i<-i
  col <- log(dataAux[[i]])
  ggp <- ggplot(dataAux, 
                aes(x = dades_netes$DEATHS, y=col)) + 
    geom_point(color = "gray30") + 
    geom_smooth(method = lm, color = "firebrick") + 
    theme_bw() + xlab("Morts") + ylab(names(dataAux)[i])
  })

}
multiplot(plotlist = histList2, cols = 4)

D’aquetes gràfiques podem veure:

  • Un increment de les variables que indiquen els danys, la intensitat i la magnitud representen un augment de morts.
  • Per les variables que fan referència a la localització (codi) i a l’any tenim una correlació negativa, pel primer cas pot ser justificable, ja que hi haurà països o zones de menor sismicitat o amb terratrèmols menys greus, però el cas dels anys no tinc explicació, (podria ser degut a la neteja de les dades que hem fet).
  • Per les altres variables no podem treure cap conclusió.

La matriu de correlacions la podem visualitzar amb la funció corrplot.

# https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html
library("corrplot")

n = c("YEAR", "DAY", "HOUR", "tsunami", "FOCAL_DEPTH", "INTENSITY", "EQ_PRIMARY", "LATITUDE", "LONGITUDE", "DEATHS", "REGION_CODE", "DAMAGE_DESCRIPTION")  

factors <- dades_netes |> 
  select(all_of(n))

# https://stackoverflow.com/questions/14674431/error-in-eigencorr-infinite-or-missing-values-in-x-when-making-a-correlat
res <- cor(factors, use = "complete.obs")
corrplot(res, 
         method = "color", 
         tl.col = "black", 
         tl.srt = 30, 
         order = "AOE",
         number.cex = 0.75,
         sig.level = 0.01, 
         addCoef.col = "black")

No veiem cap correlació significativa, a part de les comentades anteriorment, que en pugui permetre eliminar cap de les variables.

6 Construcció del conjunt de dades

No hem trobat cap correlació significativa de manera que poguem eliminar variables, per tant considerem que el nostre model no té informació que el pugui debilitar. Hem afegit una variable numèrica nova tsunami que pren valors 0 /1, i és simplement la codificació de la variable FLAG_TSUNAMI que pren valors “Yes” / “No”. Les altres variables són exactament les mateixes que hem mostrat després de la neteja.

6.1 Discretització

Abans hem vist la relació entre magnitud i profunditat, el que farem ara és discretitzar la variable FOCAL_DEPTH, la profunditat. Ho farem seguint la divisió que fa USGS dels terratrèmols:

  • Shallow (superficials): de 0 a 70 km.
  • Intermediate (intermedis): de 70 km a 300 km.
  • Deep (profunds): de 300 km a 700 km.
summary(dades_netes[,"FOCAL_DEPTH"])
##   FOCAL_DEPTH    
##  Min.   :  0.00  
##  1st Qu.: 10.00  
##  Median : 22.00  
##  Mean   : 31.67  
##  3rd Qu.: 33.00  
##  Max.   :651.00
dataAux["depth"] <- cut(dataAux$FOCAL_DEPTH, breaks = c(0, 70, 300, 700), labels = c("Shallow", "Intermediate", "Deep"))
head(dataAux$depth)
## [1] Shallow Shallow Shallow Shallow Shallow Shallow
## Levels: Shallow Intermediate Deep
plot(dataAux$depth,
     main = "Nombre de terratrèmols per segment de profunditat", 
     xlab = "Segment de profunditat", ylab = "Quantitat", 
     col = "red4")

6.2 Normalització

Tot seguit normalitzem (per la diferència), les variables per tal que cada variable contribueixi per igual al nostre anàlisi.

# Definim la funció de normalització
 nor <- function(x) {(x -min(x))/(max(x)-min(x))}

# Guardem un nou dataset normalitzat
dades_netes$type<- NULL

n = c("YEAR", "DAY", "HOUR", "tsunami", "FOCAL_DEPTH", "EQ_PRIMARY", "LATITUDE", "LONGITUDE", "DEATHS", "REGION_CODE", "DAMAGE_DESCRIPTION")

dades_netes <- dades_netes |>  
  select(all_of(n))

dades_netes_nor <- as.data.frame(lapply(dades_netes, nor))
 head(dades_netes_nor)
##          YEAR       DAY       HOUR tsunami FOCAL_DEPTH EQ_PRIMARY  LATITUDE
## 1 0.000000000 0.9333333 0.39130435       1  0.03840246  0.8513514 0.5651338
## 2 0.008333333 0.2666667 0.78260870       1  0.05069124  0.8243243 0.8224871
## 3 0.016666667 0.4000000 0.39130435       0  0.02304147  0.6486486 0.8233565
## 4 0.016666667 0.6000000 0.08695652       1  0.05069124  0.7297297 0.5912170
## 5 0.016666667 0.7000000 0.13043478       0  0.04608295  0.7567568 0.8162272
## 6 0.016666667 0.5000000 0.21739130       0  0.01382488  0.5810811 0.8242260
##   LONGITUDE       DEATHS REGION_CODE DAMAGE_DESCRIPTION
## 1 0.3148344 7.594961e-05      0.5000          0.6666667
## 2 0.8990557 5.379764e-05      0.1250          0.3333333
## 3 0.6362543 2.689882e-04      0.1875          1.0000000
## 4 0.2447166 6.325969e-03      0.5625          0.6666667
## 5 0.7136643 7.908253e-03      0.1875          1.0000000
## 6 0.7027259 1.543992e-02      0.1875          1.0000000

7 Procés de PCA

Nota: Aquest estudi el fem seguint la PAC1 i també STDA-PCA

Volem treballar amb noves característiques anomenades components principals, independents entre sí. Això ens permet representar el joc de dades en un nou sistema de coordenades, format per aquestes components principals. Aquest sistema està més ben adaptat a la distribució del joc de dades, de manera que recull millor la seva variabilitat, són noves variables que s’obtenen de combinacions lineals o barreges de les variables inicials, de manera que aquestes components principals, en particular les primeres retenen la major part de la informació.

Apliquem PCA a les columnes del nostre dataset, ho fem primer executant la funció prcomp().

pca.earthq <- prcomp(dades_netes_nor)
summary(pca.earthq)
## Importance of components:
##                           PC1    PC2    PC3    PC4    PC5     PC6     PC7
## Standard deviation     0.4314 0.3715 0.3336 0.3093 0.2845 0.22201 0.19025
## Proportion of Variance 0.2511 0.1862 0.1502 0.1290 0.1092 0.06649 0.04883
## Cumulative Proportion  0.2511 0.4373 0.5874 0.7165 0.8256 0.89214 0.94097
##                            PC8     PC9    PC10    PC11
## Standard deviation     0.16748 0.10051 0.05801 0.04736
## Proportion of Variance 0.03784 0.01363 0.00454 0.00303
## Cumulative Proportion  0.97881 0.99243 0.99697 1.00000

La funció summary(), ens mostra la proporció de variància aplicada al conjunt total de cada atribut. Podem veure que l’atribut 1 explica el 0.2511 de variabilitat del total de dades, en canvi l’atribut 7 explica només el 0.04883.

El següent histograma mostra el pes de cada atribut sobre el conjunt total de dades:

library('factoextra')

#Els valors propis corresponen a la quantitat de variació explicada per cada component principal (PC).
ev= get_eig(pca.earthq)
ev
##         eigenvalue variance.percent cumulative.variance.percent
## Dim.1  0.186125657       25.1090145                    25.10901
## Dim.2  0.138024978       18.6200615                    43.72908
## Dim.3  0.111304141       15.0153253                    58.74440
## Dim.4  0.095655537       12.9042728                    71.64867
## Dim.5  0.080915617       10.9158051                    82.56448
## Dim.6  0.049290469        6.6494599                    89.21394
## Dim.7  0.036194481        4.8827645                    94.09670
## Dim.8  0.028049490        3.7839761                    97.88068
## Dim.9  0.010101961        1.3627905                    99.24347
## Dim.10 0.003364748        0.4539165                    99.69739
## Dim.11 0.002243182        0.3026132                   100.00000
fviz_eig(pca.earthq)

Seguint el guió de la PAC1, volem utilitzar el mètode de Kàiser per decidir quina de les variables obtingudes farem servir. Aquest criteri manté totes les variables la variància de les quals és superior a 1.

# Calculem la variància dels components principals a partir de la desviació estàndard

var_earthq <- pca.earthq$sdev^2

var_earthq
##  [1] 0.186125657 0.138024978 0.111304141 0.095655537 0.080915617 0.049290469
##  [7] 0.036194481 0.028049490 0.010101961 0.003364748 0.002243182

Aquests valors obtinguts no ens deixen clar quines components hem de considerar. Per ajudar-nos a triar-les, escalem les dades, cosa que no havíem fet abans i tornem a calcular la variància.

# Escalem les dades
earthq_scale <- scale(dades_netes_nor)
# Calculem les components principals
pca.earthq_scale <- prcomp(earthq_scale)
# Mostramos la varianza de dichas variables:
var_earthq_scale <- pca.earthq_scale$sdev^2
head(var_earthq_scale)
## [1] 2.0645709 1.6218845 1.1125793 1.0508882 1.0385397 0.9398405

Després d’analitzar la variància i aplicant el criteri de Kàiser ens quedarem amb les components principals 1, 2, 3, 4 i 5 que són les superiors a 1. Aquest criteri té el problema de sobreestimar el nombre de factors, però malgrat això és el que aplicarem per analitzar els resultats.

Mostrem l’histograma de percentatge de variància explicat amb les dades escalades:

fviz_eig(pca.earthq_scale)

ev = get_eig(pca.earthq_scale)
ev
##        eigenvalue variance.percent cumulative.variance.percent
## Dim.1   2.0645709        18.768827                    18.76883
## Dim.2   1.6218845        14.744405                    33.51323
## Dim.3   1.1125793        10.114357                    43.62759
## Dim.4   1.0508882         9.553529                    53.18112
## Dim.5   1.0385397         9.441270                    62.62239
## Dim.6   0.9398405         8.544005                    71.16639
## Dim.7   0.8861009         8.055463                    79.22185
## Dim.8   0.7756779         7.051617                    86.27347
## Dim.9   0.5961553         5.419594                    91.69307
## Dim.10  0.5322558         4.838689                    96.53175
## Dim.11  0.3815070         3.468245                   100.00000

Aquesta taula s’interpreta de la següent manera, la suma de tots els valors propis és 11. Així considerant el primer valor propi, 2.0645709 i dividint-lo per 11, obtenim el 18.768827 percent de la tercera columna i la darrera és la suma dels valors anteriors acumulats.

Aquests valors propis es poden utilitzar per determinar el nombre de components principals a retenir després de la PCA (Kaiser 1961):

Observem que les 5 primeres components principals expliquen el 62.6% de la variació. És un valor una mica baix, però en aquest cas, per les dades que tenim ho donem per bo.

Així, després d’aplicar el mètode Kàiser ens quedem amb 5 components principals. Per obtenir els resultats d’aquestes 5 variables del que ens ha donat la sortida del PCA, fem servir la funció get_pca_var() (factoextra package). Aquesta funció ens dona una llista de matrius que contenen el resultat de les variables considerades (coordenades, correlació entre variables i eixos, cosinus al quadrat i contribucions).

var <- get_pca_var(pca.earthq_scale)
var
## Principal Component Analysis Results for variables
##  ===================================================
##   Name       Description                                    
## 1 "$coord"   "Coordinates for the variables"                
## 2 "$cor"     "Correlations between variables and dimensions"
## 3 "$cos2"    "Cos2 for the variables"                       
## 4 "$contrib" "contributions of the variables"

Els components de get_pca_var() es poden utilitzar en el diagrama de variables de la següent manera:

#Utilitzem els 5 components principals trobats abans
head(var$coord[,1:5],11)
##                          Dim.1       Dim.2       Dim.3       Dim.4       Dim.5
## YEAR               -0.49726751 -0.26462329  0.22745564  0.25608947 -0.34108621
## DAY                 0.07749664  0.02575068 -0.20154527 -0.67088882 -0.26413864
## HOUR                0.04217844  0.03367169  0.05360339 -0.07727443 -0.86173371
## tsunami             0.67205414  0.22971898  0.24013393  0.24746499 -0.01419688
## FOCAL_DEPTH         0.31095625 -0.16922316  0.40773756 -0.61222775  0.12059281
## EQ_PRIMARY          0.78900900  0.26567688  0.25134400  0.02580474  0.01789808
## LATITUDE           -0.42314110  0.53895037 -0.35378868 -0.21624743  0.18312465
## LONGITUDE          -0.36993791  0.52587765  0.53467986  0.03843890  0.01096036
## DEATHS              0.19837872  0.40833002 -0.18684736  0.14628881 -0.23133026
## REGION_CODE         0.44251807 -0.60327882 -0.35527994  0.13439894 -0.00372567
## DAMAGE_DESCRIPTION  0.29588280  0.54806669 -0.39758852  0.06936600 -0.08717138

7.1 Qualitat de la representació

La qualitat de representació de les variables en el mapa de factors s’anomena cos2 (cosinus quadrat, coordenades quadrades). Podem accedir al cos2 de la següent manera:

head(var$cos2[,1:5],11)
##                          Dim.1        Dim.2       Dim.3        Dim.4
## YEAR               0.247274973 0.0700254873 0.051736069 0.0655818182
## DAY                0.006005729 0.0006630974 0.040620496 0.4500918103
## HOUR               0.001779021 0.0011337825 0.002873324 0.0059713377
## tsunami            0.451656771 0.0527708094 0.057664303 0.0612389197
## FOCAL_DEPTH        0.096693789 0.0286364784 0.166249918 0.3748228133
## EQ_PRIMARY         0.622535208 0.0705842045 0.063173807 0.0006658845
## LATITUDE           0.179048389 0.2904675000 0.125166433 0.0467629512
## LONGITUDE          0.136854054 0.2765473012 0.285882554 0.0014775493
## DEATHS             0.039354115 0.1667334030 0.034911937 0.0214004159
## REGION_CODE        0.195822246 0.3639453363 0.126223834 0.0180630749
## DAMAGE_DESCRIPTION 0.087546630 0.3003771000 0.158076628 0.0048116422
##                           Dim.5
## YEAR               1.163398e-01
## DAY                6.976922e-02
## HOUR               7.425850e-01
## tsunami            2.015515e-04
## FOCAL_DEPTH        1.454263e-02
## EQ_PRIMARY         3.203412e-04
## LATITUDE           3.353464e-02
## LONGITUDE          1.201296e-04
## DEATHS             5.351369e-02
## REGION_CODE        1.388061e-05
## DAMAGE_DESCRIPTION 7.598849e-03

Podem visualitzar el cos2 de les variables en totes les dimensions amb el paquet corrplot:

corrplot(var$cos2[,1:5], is.corr=FALSE)

També podem crear un diagrama de barres de variables cos2 amb la funció fviz_cos2():

fviz_cos2(pca.earthq_scale, choice = "var", axes = 1:2)

Nota:

  • Un cos2 elevat indica una bona representació de la variable en el component principal. En aquest cas, la variable es dibuixada prop de la circumferència del cercle de correlació.

  • Un cos2 baix indica que la variable no està perfectament representada pels PC. En aquest cas, la variable està a prop del centre del cercle.

Per a una variable donada, la suma del cos2 de tots els components principals és igual a 1.

Si una variable està perfectament representada per només dos components principals (Dim.1 i Dim.2), la suma del cos2 en aquests dos PCs és igual a un. En aquest cas les variables es col·locaran en el cercle de correlacions.

Per a algunes de les variables, poden ser necessaris més de 2 components per representar perfectament les dades. En aquest cas les variables es situen dins del cercle de correlacions.

Resumint:

  • Els valors de cos2 s’utilitzen per estimar la qualitat de la representació
  • Com més propera estigui una variable al cercle de correlacions, millor serà la seva representació al mapa de factors (i més important és interpretar aquests components)
  • Les variables que estan properes al centre de la trama són menys importants per als primers components.
fviz_pca_var(pca.earthq_scale,
             col.var = "cos2", 
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE     
             )

7.2 Contribució

Les contribucions de les variables en la comptabilització de la variabilitat d’un determinat component principal s’expressen en percentatge.

  • Les variables que estan correlacionades amb PC1 (és a dir, Dim.1) i PC2 (és a dir, Dim.2) són les més importants per explicar la variabilitat en el conjunt de dades.

  • Les variables que no estan correlacionades amb cap PC o amb les darreres dimensions són variables amb una contribució baixa i es poden eliminar per simplificar l’anàlisi global.

La contribució de les variables es pot extreure de la següent manera:

head(var$contrib[,1:5],11)
##                          Dim.1       Dim.2      Dim.3       Dim.4        Dim.5
## YEAR               11.97706360  4.31753848  4.6501017  6.24060839 11.202248904
## DAY                 0.29089477  0.04088438  3.6510202 42.82965619  6.718011866
## HOUR                0.08616904  0.06990526  0.2582579  0.56821816 71.502800880
## tsunami            21.87654422  3.25367247  5.1829387  5.82734859  0.019407200
## FOCAL_DEPTH         4.68348110  1.76562994 14.9427477 35.66723912  1.400295553
## EQ_PRIMARY         30.15324882  4.35198712  5.6781397  0.06336397  0.030845348
## LATITUDE            8.67242616 17.90925926 11.2501134  4.44985018  3.229018205
## LONGITUDE           6.62869231 17.05098614 25.6954766  0.14060004  0.011567165
## DEATHS              1.90616436 10.28022667  3.1379279  2.03641220  5.152782374
## REGION_CODE         9.48488828 22.43965808 11.3451539  1.71883885  0.001336551
## DAMAGE_DESCRIPTION  4.24042734 18.52025221 14.2081223  0.45786432  0.731685955

Quant més gran sigui el valor de la contribució, més aportació de la variable hi haurà al component.

corrplot(var$contrib[,1:5], is.corr=FALSE)

Les variables més importants (que més contribueixen) es poden ressaltar a la gràfica de correlació de la següent manera:

fviz_pca_var(pca.earthq_scale, col.var = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07")
             )

Les variables que apunten al mateix costat del gràfic estan correlacionades de manera positiva, i les que apunten a costats oposats del gràfic ho estan negativament.
Així, podem fer les següents observacions:

  • La variable EQ_PRIMARY, (magnitud, és la que més aporta), està força correlacionada amb tsunami, a la matriu de correlacions s’ha vist que tenia un valor de 0.52. De fet és normal ja que aquesta darrera indica la presència d’un terratrèmol d’una certa magnitud. Ara bé, no vol dir que sempre hi ha tsunamis, això depèn d’altres factors, principalment de si és una zona marítima.
  • Una altra correlació és entre DEATHS i DAMAGE_DESCRIPTION, és a dir els danys i els morts per terratrèmols estan relacionats.
  • Les quatre variables anteriors apunten a la mateixa direcció i per tant poden formar un indicador de la gravetat i destrucció d’un terratrèmol.
  • Si ara ens fixem amb FOCAL_DEPTH i REGION_CODE, també tenen una correlació positiva. Potser podria indicar alguna característica geològica de les zones on han passat els sismes.
  • YEAR està correlacionada d’una manera negativa amb les 4 anteriors, en principi no treurem conclusions clares ja que hem fet una neteja exhaustiva, potser pot tenir relació amb que amb anys més recents tenim més dades.
  • En quant a LONGITUDE i LATITUDE, formen una parella que ens serveixen per situar els sismes en el mapa.

8 Conclusions

Hem treballat amb un dataset, Global Significant Earthquake Database from 2150BC, que és un llistat d’uns 6000 terratrèmols significatius que van del 2150 BC fins a l’actualitat. Aquests terratrèmols (registres) tenen un identificador únic i una sèrie de variables que el descriuen: variables temporals que en donen informació de quan van passar els sismes, variables que descriuen les característiques de terratrèmol, variables geogràfiques i variables que descriuen els efectes i danys sobre les persones i les infraestructures.

Al ser un dataset amb una gran presència de valors nuls, hem hagut de descartar una gran quantitat de registres, de fet hem considerat la sismicitat significativa a partir de l’any 1900, que és quan l’increment de les dades disponibles va augmentant, i hem netejat les variables que més ens han interessat, mantenint la resta. Això ha minvat la qualitat de les dades, tot i així ens ha permès seguir amb l’estudi. (Val a dir que amb paciència es podria cercar o construir un dataset més complet).
Ens hem quedat al final amb 1135 registres, (tot i que per alguna gràfica, en particular els mapes hem fet servir més punts per posar de manifest algunes estructures o regularitats).

Hem vist que esl terratrèmols estan distribuïts seguint els límits de les plaques tectòniques i que alguns d’ells, els que estan en zones marítimes desencadenen tsunamis. Aquests provoquen un augment de morts, però com que el nombre de terratrèmols sense tsunami és molt més superior, la mortalitat també ho és.
S’ha vist que els terratrèmols amb tsunami generalment tenen magnitud superior a 6 aproximadament.

També hem vist que l’increment de la magnitud del terratrèmol incrementa el nombre de víctimes, però a la vegada s’observa la presència de magnituds elevades amb pocs morts, això ens indica que la densitat de població pot ser un factor important.

Visualitzant la relació entre la magnitud i la profunditat hem pogut determinar que els terratrèmols a una profunditat superior a 300 km, la magnitud és més gran que 7.

Posant noms als països amb més nombre de terratrèmols i les zones amb els més forts, mortals i destructius hem deduït que factors com la densitat de població i riquesa i desenvolupament dels països afectats poden influir en el nombre de morts i danys ocasionats.

Finalment, amb l’estudi PCA hem vist una sèrie de relacions que hem comentat i també que amb una nova variable combinació d’altres amb una correlació inicial, es podria considerar com un índex de devastació o destrucció d’un terratrèmol (significatiu).